Calculation of the incremental stress-strain relation of a polygonal packing 
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The constitutive relation of the quasi-static deformation on two dimensional packed samples 
of polygons is calculated using molecular dynamic simulations. The stress values at which the 
system remains stable are bounded by a failure surface, that shows a power law dependence on 
the pressure. Below the failure surface, non linear elasticity and plastic deformation are obtained, 
which are evaluated in the framework of the incremental linear theory. The results shows that the 
stiffness tensor can be directly related to the micro-contact rearrangements. The plasticity obeys a 
non-associated flow rule, with a plastic limit surface that does not agree with the failure surface. 
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I. INTRODUCTION 

The non linear and irreversible behavior of soils has 
been described by different constitutive theories [|l],D. 
Here the stress - strain relation is postulated using a cer- 
tain number of material parameters that are measured 
in experimental tests. In practice, the constitutive re- 
lations can be constructed directly, taking samples with 
the same macroscopical state, and measuring in each one 
the incremental strain that results from the application 
of a specific stress increment However, such tests are 
difficult to perform, because they require the fabrication 
of many samples with identical material properties. 

Numerical simulations result as an alternative to the 
solution of this problem. They allow to develop dif- 
ferent tests on identically generated samples and pro- 
vide detailed information about micro-mechanical rear- 
rangements during the loading process. Usually, disks or 
spheres are used to capture the granularity of the sample 
Although the simplicity of their geometry allows 
to reduce the computer time of calculations, they do not 
provide a detailed description of realistic granular tex- 
tures. 

We present here a two dimensional discrete model that 
takes into account the diversity of shapes of the grains 
in the soils. The granular samples consist of randomly 
generated polygons. As presented by Tillemans The 
interaction between the polygons could be handled let- 
ting the polygons interpenetrate each other and calcu- 
lating the force as a function of their overlap. This 
approach has been successfully applied to model differ- 
ent processes, like fragmentation and strain localiza- 
tion p[. A suitable contact force law is introduced in 
Sect 



[I A 



that attempts to combine the Hertz contact 
law with the Coulomb friction criterion for quasi-static 
deformation. 

The incremental stress-strain relation was calculated 
performing different stress increments on the same sam- 
ple, and measuring the corresponding strain response. 
The stress is applied on the boundary through a flexible 



membrane that surrounds the sample. The modeling of 
such a membrane, whose details are given in Sect. II B, 
results more complex than rigid walls. However, it re- 
sults more advantageous than walls, because it allows to 
implement a stress-controlled condition without any re- 
striction in the deformation of the boundary. 



II. MODEL 

The polygons of this model are generated using a sim- 
ple version of the Voronoi tessellation: First we set a ran- 
dom point in each cell of a regular square lattice. Then 
each polygon is constructed assigning to each point that 
part of the plane that is nearer to it than to any other 
point. Each polygon is subjected to interparticle contact 
forces and boundary forces that are inserted in Newton's 
equation of motion. 



A. Contact force 

Usually, the interaction between two solid bodies in 
contact is described by a force applied on the flattened 
contact surface between them. Given two polygons in 
contact, such surface is obtained from the geometrical 
construction shown in Fig. |l|. The points Ci and C2 re- 
sult from the intersection between the edges of the poly- 
gons. The contact surface is taken as the segment that 
lies between those points. The vector S = C1C2 defines 
an intrinsic coordinate system at the contact (t, h) , where 
t = >5'/|S'| and n is perpendicular to it. The deformation 
length is given by S — a/\S\ where a is the overlap area 
between the polygons. £ is the branch vector, that con- 
nects the center of mass of the polygon to the point of 
application of the contact force, that is supposed to be 
the center of mass of the overlap area. 

The normal elastic force is taken proportional to the 
deformation length as = fc„J; the tangential force 
is calculated from the simplified Coulomb friction law, 
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with a single friction coefficient = l^d = ^J^■ Here 
Us is the static and /x^ the dynamic friction coeflicient. 
This tangential force is implemented by an elastic spring 
f^ = —kt^, where ^ grows linearly with the tangential 
displacement of the contact, whenever \f^ \ < /i/^. We 
used the straightforward calculation of ^ proposed by 
Brendel [|: 

f vn{t')Q{mt')-ti\ft{t')\)dt' (1) 

Jo 

where @ is the Heaviside function and v is the relative 
velocity at the contact, that depends on the linear veloc- 
ity Vi and angular velocity uji of the particles in contact 
according to: 

V = Vi — Vj — iSi X ii + LOj X ij (2) 
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FIG. 1. Contact surface as defined from the geometry of 
overlap. 

B. Boundary Forces 

Let us now discuss how to apply the stress on the sam- 
ple. One way to do that, would be to apply a perpendic- 
ular force on each edge of the polygons belonging to the 
external contour of the sample. Actually this does not 
work, because this force will act on all the fjords of the 
boundary. It produces an uncontrollable growth of cracks 
that with time ends up destroying the sample. Thus, re- 
sults necessary to introduce a flexible membrane in order 
to restrict the boundary points that are subjected to the 
external stress. 

The algorithm to identify the boundary is rather sim- 
ple. The lowest vertex -p from all the polygons of the 
sample is chosen as the first point of the boundary list 
&i. In Fig. |2| P is the polygon that contains p, and 
q G P n Q is the first intersection point between the 
polygons P and Q in counterclockwise orientation with 
respect to p. Starting from p, the vertices of P in coun- 
terclockwise orientation are included in the boundary list 



until q is reached. Next, q is included in the boundary 
list. Then, the vertices of Q between q and the next 
intersection point r G Q fl P in the counterclockwise ori- 
entation are included into the list. The same procedure 
is applied until one reaches the lowest vertex p again. 
This is a very fast algorithm, because it only makes use 
of the contact points between the polygons, which are 
previously calculated to obtain the contact force. 
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FIG. 2. Algorithm used to find the boundary. 

The set of points that are in contact with the mem- 
brane are selected using a recursive algorithm. It is ini- 
tialized with the vertices of the smallest convex polygon 
that encloses the boundary, (see Fig. ||). The lowest 
point of the boundary is selected as the first vertex of the 
polygon mi h\. The second one mi is the boundary 
point hi that minimizes the angle L{h\hi) with respect to 
the horizontal. The third one TO3 is the boundary point 
hi such that the angle /.(mibi, mi mi) is minimal. The 
algorithm is recursively applied until the lowest vertex 
TOi is reached again. 

The points of the boundary are iteratively included 
in the list m^ using the bending criterion proposed by 
Astr0m [0 : For each pair of consecutive vertices of the 
membrane rrii = hi and rrii^i = hj we choose that point 
from the subset {hk}i<k<j that maximizes the bending 
angle 9h = /_{hkhi,hkhj). This point is included into the 
list, whenever 6b > dth- Here 9th is a threshold angle for 
bending. This algorithm is repeatedly applied until there 
are not more points satisfying such bending condition. 

The final result gives a set of segments {mirdi^i} ly- 
ing on the boundary of the sample. In order to apply 
the boundary forces, those segments are divided into two 
groups: A- type segments are those that coincide with an 
edge of a boundary polygon; B-type segments connect 
the vertices of two different boundary polygons. 
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FIG. 3. Membrane obtained with threshold bending angle 
8th ~ 1", 37r/4, 7r/2 and 7r/4, the first one corresponds to the 
minimum convex polygon that encloses the sample. 

On each segment of the membrane ff^rjl^_^_i a force 
fi = (JiNi is appUed, where Ui is the local stress and 
Ni is the 90^^ counterclockwise rotation of fr^frlij^i. This 
force is transmitted to the polygons in contact with it: 
if the segment is A-type, this force is applied in its mid- 
point; If the segment is B-type, half of the force is applied 
at each one of the vertices connected by this segment. 



C. Molecular dynamic simulation 

Before we implement the numerical solution of New- 
ton's equations it is convenient to make a dimensional 
analysis of the parameters. In such way we can keep the 
scale invariance of the model and reduce the parameter to 
a minimum of adimensional constants. All the polygons 
are supposed to have the same density. The mass rrii of 
each polygon is measured in units of the mean mass toq 
of the Voronoi tessellation. The time is measured in frac- 
tions of the total loading time t^. The evolution of the 
position Xi and the orientation ipi of the i — th polygon 
is governed by the equations of motion: 



A rriiXi -f 



A^T,^, ^4 X /f + ^ -^^r X /f = 

c Cb 



(3) 



The sums go over all those particles and boundary seg- 
ments that are in contact with the i — th polygon. The 
interparticle contact forces and boundary forces are 
given by: 



I, 



b _ 



+ m - A7m<)t? 



Here and are the deformation length and the tan- 
gential displacement of the contact, which were defined 
in Sect. |II A|; cr? is the stress applied on the boundary 

Artificial viscous terms 



II B 



II A| ; 

segment T^, defined in Sect, 
must be included in Eq. (^) to keep the stability of the 
numerical solution and reduce the acoustic waves gener- 
ated during the loading process, v'^ is the relative velocity 
at the contact (Eq. (Q)) and m = (l/m^ -I- l/mj)~^ is 
the effective mass of the two polygons in contact. 

There are four microscopic parameters in the model: 
the viscosity 7, the ratio A = tg/to between the character- 
istic period of oscillation ts = i/fcn/mo and the loading 
time to, the friction coefficient and the ratio C = kt/kn 
between the tangential kt and normal fc„ stiffness of the 
interparticle contacts. 

The viscosity factor 7 is related to the normal restitu- 
tion coefficient |9). It was taken large enough to have a 
high dissipation, but not too large to keep the numeri- 
cal stability of the method. The ratio A was chosen small 
enough in order to to avoid rate-dependence in the strain 
response, as corresponds to the quasi-static approxima- 
tion. Technically, that is done by looking for the value of 
A such that a reduction of it by half makes a change of 
the strain response less than 5%. 

The two parameters ^ and fj, determine the constitu- 
tive response of the system. For example, the micro- 
mechanical analysis of the strain response shows that the 
Young's modulus and Poisson's ratio depend on ( [0. 
In the other hand, /i can be directly related to the fric- 
tion angle of the material |jl^ . Although the study of the 
dependence of the constitutive response on those param- 
eters is an important point, such quantities have been 
kept fixed in this work. 

The boundary conditions yield more dimensional pa- 
rameters. The initial hight Hq and width Wq of the sam- 
ple, and the characteristic length £0 of the polygons define 
two geometrical parameters, which are the shape ratio 
Wq/Hq and the granularity £o/Ho of the sample. 

In order to keep overlaps much smaller than the char- 
acteristic area of the polygons. The ratio Ui/kn between 
the stress applied on the membrane and the stiffness of 
the contacts is restricted to small values. This was imple- 
mented by fixing the contact stiffness to a value closed to 
the experimental granular stiffness /c„ = l&OMPa. Then 
the stress is chosen in such way that it does not excess 
1% of this value. 



(4) 



adimensional variable 


ratio 


default value 


viscosity 


7 


0.1 


friction coefficient 




0.25 


time's ratio 


A = ts/to 


8.0 X 10-4 


stiffness ratio 


^ = kt/kn 


0.33 


granularity 


io/Ho 


0.1 


shape ratio 


Wo/Ho 


1.0 


bending angle 




0.257r 
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III. STRESS-STRAIN CALCULATION 

A. Theorical background 

The macroscopic state of the system is characterized 
by the stress tensor and the void radio e. The area frac- 
tion of voids in the sample defines the void ratio. Initially 
Co = due to the Voronoi tessellation used. The stress 
controlled test was restricted to stress states without-off 
diagonal components. The diagonal component, the ax- 
ial Si and lateral S3 stress, define the stress vector: 



(5) 



where p and q are the pressure and the shear stress. The 
domain of admissible stresses is bounded by the failure 
surface. When the system reaches this surface it becomes 
unstable and fails. 

Before failure, the constitutive behavior can be ob- 
tained performing small changes in the stress and evalu- 
ating the resultant deformation. An infinitesimal change 
of the stress vector da produces an infinitesimal defor- 
mation of the sample, which is given by a change of 
height dH and width dW. This defines the axial strain 
dei ~ dH/H and lateral strain des — dW/W increments. 
The volumetric strain dey and the shear strain de^ incre- 
ments define the incremental strain vector: 



p 


1 


■<5i 


+ S3 ' 


_ q _ 


~ 2 


.'^1 


-S3_ 



dl- 



dey 




dei + des 


de~f 




dei — des 



(6) 



Each state of the sample is related to a single point 
in the stress space, and the quasi-static evolution of the 
system is represented by the movement of this point in 
the stress space. The constitutive relation is formulated 
taking the incremental strain as a function of the incre- 
mental stress and the stress state. 



de — J- {da, a) 



(7) 



If there is no rate - dependence in the constitutive 
equation, !F{da) is an homogeneous function of degree 
one. In this case, the application of the Euler identity ||] 
shows that Eq. (|^) can be reduced to. 



de = M{§, a)da 



(8) 



Where 9 is the unitary vector defining a specific direction 
in the stress space: 



da 



cost) 
sin 9 



\da\ = \/ d-p^ + dcp- 



(9) 



The constitutive relation results from the calculation 
of de{9), where each value of 9 is related to a particular 
mode of loading. Some special modes are listed in the 
table: 



0° 


isotropic compression 


dp > dq = 


45° 


axial loading 


dai > da^ = 


90° 


pure shear 


dp = dq > 


135° 


lateral loading 


dai = das > 


180° 


isotropic expansion 


dp <0 dq = 


225° 


axial stretching 


dai < das = 


270° 


pure shear 


dp^ dq <0 


315° 


lateral stretching 


dai = das < 



The relation (||) has been proposed by Darve , and 
it contains all the possible constitutive equations. In or- 
der to interpret our particular results, it is convenient 
to make some approximations: First, if the load incre- 
ments are taken small enough, the tensor M{9) can be 
supposed to be lineal in each stress direction. Then, we 
assume that the strain can be separated in an elastic (re- 
coverable) and a plastic (unrecoverable) component: 



de = de'' + deP 



de" = D{a)da 



deP = J{9,d)da 



(10) 



(11) 



(12) 



Here, D ^ defines the stiffness tensor, and J = M — D 
the flow rule of plasticity, which results from the calcu- 
lation of ^6^(6*) and deP{e). 




FIG. 4. Axial stress tJi — p + q and lateral stress as = p~q 
in a stress controlled test. They are applied on the boundary 
of the tessellated sample of polygons. 
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B. The method 

A numerical method has been implemented in order 
to find the elastic de'^ and plastic de^ components of the 
strain as function of the stress state a and the stress di- 
rection 9. Fig. H shows the three steps of the procedure: 

1) The sample is driven to the stress state a. First, 
it is isotropically compressed until it reaches the stress 
value ^1 — S3 — p — q. Next, it is subjected to axial 
loading, in order to increase the axial stress Si until p + q 
(see Fig. ^. When the stress state a = [p q]'^ is reached, 
{A^ being the transpose of A) the sample is allowed to 
relax. 

2) Loading the sample from a to a + da the strain in- 
crement de is obtained. This procedure is implemented 
choosing different stress directions according to Eq. (jp). 
Here the stress modulus is fixed to \dd-\ — 10~'^p. 

3) The sample is unloaded until the original stress state 
a is reached. Then one finds a remaining strain de^, that 
corresponds to the plastic component of the incremen- 
tal strain. Since the stress increments are taken small 
enough, the unloaded stress-strain path is practically 
elastic. Thus, the difference de'^ = de — de^ represents 
the elastic component of the strain. 




dp 



FIG. 5. Procedure to obtain the constitutive behavior: 1) 
The sample is driven to the stress state a, with pressure p 
and shear stress q. 2) It is loaded from ct to ct -I- dtr 3) It is 
unloaded to the original stress state a. 

One could be concerned about the dependence of the 
strain response on the way how the stress state is reached. 
We found that there is not remarkable dependence of 
the strain response on the stress path, whenever the 
stress components are quasi - static and monotonically 
increased. Otherwise, a strong reduction in the plas- 
tic component of the strain is observed. In fact, when 
the plastic response is calculated after the sample is un- 
loaded, the plasticity results smaller than that one cal- 
culated after a monotonic load. Furthermore, there is 



no plastic component in the strain response when elas- 
tic waves are previously generated in the sample. Those 
memory effects suggest that the plastic component of the 
strain depend on the history of the deformation, and is 
kept unchanged only if the sample is subjected to quasi- 
static and monotonic load. 

Fig. H shows the load-unload paths and the corre- 
sponding strain response. They were taken from a stress- 
state with q = Q.5p. The end of the load paths in the 
stress space map into a strain envelope response d€{9) in 
the strain space. Likewise, the end of the unload paths 
map into a plastic envelope response deP{9). The yield 
direction (j) can be found from this response, as the di- 
rection in the stress space where the plastic response is 
maximal. The flow rule can be obtained taking the di- 
rection ip of the maximal plastic response in the strain 
space. These angles do not agree, that reveals the ne- 
cessity to analyze this behavior in the framework of the 
non-associated theory of plasticity (see Sect. IV C). 



, X 10 




to 
T3 







X 10"' 




dp/p 


1 


135°, 


-.-90 














180° \ 


















- 0° 


< 

225° / 


i \ ■> 

270° 


315° 



de 



1 2 

X 10'' 



FIG. 6. Stress - strain relation resulting from the load - 
unload test. Dotted lines represent the paths in the stress 
and strain spaces. The dash-dot line gives the strain envelope 
response and the solid line is the plastic envelope response. 
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IV. CONSTITUTIVE RELATION 

Fig.0 summarizes the global elasto-plastic behavior. 
The elastic response, calculated from Eq. (|o|), has a cen- 
tered ellipse as envelope response. This can be related to 
the micro-contact structure using a local li near r elation 
in each point of the stress space (see Sect. IV b| ). The 
solid line represents the failure surface, which separates 
the stable states of the unstable ones (see Sect. IV A ). 
The plastic envelope response is almost on a straight line. 
The modulus and the orientation of this envelope depend 
on the stress state through a certain number of material 
parameters, which are given in Sect. [VC. All the quan- 
tities obtained in this section have been calculated from 
the average over five different samples of 10 x 10 particles 
each one. 
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FIG. 7. Elastic response de"^ and plastic response de^ re- 
sulting from the application of different loading modes with 
\da\ — 10~*p. The solid line represents the failure surface. 



A. Failure surface 

The failure line was calculated looking for the values of 
stress for which the system becomes unstable: for each 
pressure p, there is a critical shear stress qdp), below 
which the sample reaches a stable state with an expo- 
nential decay of its kinetic energy. For shear stress values 
above the critical one, the sample develops an instability 
and fails. Fig. || shows the interface between these two 
stress states, that can be accurately fitted by the power 
law: 



qc_ 

Po 



/3 



(13) 



Here po = l.OMPa is the reference pressure, and fi* 
0.78 ± 0.03 is the Mohr-Coulomb friction coefficient 



The power law dependence on the pressure, with expo- 
nent (3 = 0.92 ± 0.02 implies a significant deviation from 
the Mohr-Coulomb theory. Moreover, the empirical cri- 
teria of failure for most rocks |l5[ shows a power law 
dependence of the form of Eq. (|13|). It seems that ad- 
ditional features beyond the Mohr-Coulomb analysis are 
taking place when the sample fails, that will be discussed 



in sect. IV C 
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FIG. 8. Failure surface. The continuous line represents the 
power law fit. 



B. stiffness 

Hooke's law of elasticity states that the stiffness tensor 
of isotropic materials can be written in terms of two ma- 
terial parameters: the Young's modulus E and the Pois- 
son's ratio v. However, the isotropy is not fulfilled when 
the stress state is far from the hydrostatic axis. Indeed, 
numerical simulations | p^ , pTt and photo-elastic experi- 
ments on granular materials show that the loading 
induces a significant departure from isotropy in the con- 
tact network. 

The anisotropy of the granular sample can be charac- 
terized by the distribution of the micro-contact normal 
vectors hi (see Fig. Our numerical simulations show 
that the structural changes of micro-contacts are princi- 
pally due to the opening of contacts whose normal vectors 
are nearly aligned around the direction perpendicular to 
the load. Let us call N(Lp)^(p the number of contacts per 
particle, oriented between the angles and (p-t- A(p, mea- 
sured with respect to the direction in which the sample is 
loaded. The lowest order of anisotropy can be described 
by: 

iV((p) = -^[7V-|-(iVo-7V)cos(2<^)] (14) 

Here iV is the average coordination number of the poly- 
gons, whose initial value iVo = 6.0 reduces as the load is 
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increased. Fig. || shows this reduction. A critical hne is 
found around q = 0.12p, below which there are no struc- 
tural changes in the contact network. Above this limit an 
induced anisotropy arises due to opened contacts whose 
amount follows a power law dependence. 

In order to describe the effect of the anisotropy in the 
elastic response we proceed as follows: first, an additional 
parameter a is included in Hooke's law 



" del ' 
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del_ 





1 



- a 
V 1 





da I 




da^ 



(15) 



Then, these three parameters are supposed to be de- 
pended on the internal damage parameter d: 



No~N 

Nn 



(16) 



The tensor D defined in Eq. (|T^) is calculated from Eq. 
( |l5| ) using the definition of the stress and strain vectors 
given in Eqs. (||) and (^). One obtains: 



D = 



l-v - 
—a 1 



(17) 



The diagonal components of this tensor are respectively 
the inverse of the bulk modulus and of the shear modulus. 
The non-diagonal component results from the anisotropy 
of the sample, and it couples the compression mode with 
the shearing deformation. These three variables are cal- 
culated from the elastic response di" {0) by the introduc- 
tion of the following function: 



R{e) 



da^de" 

\da\^ 



(18) 



by substitution of Eqs. ( |ll| ) and (g) into Eq. (18), one 
sees that R is the quadratic form of D: 

R{0) = e^^De cos(26l) - a sin(26l)] (19) 

Using this equation, the components of D can be evalu- 
ated as the Fourier coefficients of R: 



1 



Ri^e)de 



a 
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4vr JO 
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E 
'2^ 



2ir 



R{e) cos{2e)d0 



R{e) sin{2e)de 



(20) 



(21) 



(22) 



Figs. |T^, and ^ show the results of the calculation 
of the Young's modulus E, the Poisson's ratio v and the 
anisotropy factor a. Below the limit of isotropy, Hooke's 
law can be applied: E Eq, u k, and a « 0. On the 
other hand, above the limit of isotropy a reduction of the 
Young's modulus is found, along with an increase of the 
Poisson's ratio and the anisotropy factor. The functional 



dependence of those parameters with the internal dam- 
age parameter d are evaluated developing their Taylor's 
series around d = 0: 



E{d) = E{0) 
a(d) ^ a(0) - 
v{d) = - 



-E'{0)d 
-a'{0)d- 
iy'{Q)d^ 



f O {d') 

-o{d^) 

iy"iO)d^ 



(23) 



0{d^ 



The coefficients of this expansion are calculated from 
the best fitting of those expansions. Figs. |l^ and |l^ 
show that the linear approximation is good enough to 
reproduce the Young's modulus and the anisotropy fac- 
tor. The fit of the Poisson's ratio, however, requires the 
inclusion of a quadratic approximation, implying that it 
has a non-linear dependence on the damage parameter 
(Fig. |l|). 



C. Plastic Flow 



The formulation of the non-associated theory of plas- 
ticity requires the evaluation of three material functions: 
the yield direction 0, the flow direction ip, and the plastic 
modulus h. These quantities can be calculated from the 
plastic response deP{9), as follows: 

The yield direction is given by the incremental stress 
direction (p with maximal plastic response: 



\dlP{(l))\ =max|rf?(6l)| 



(24) 



The flow direction is defined from the orientation of the 
plastic response at its maximum value: 



deP 

i(j = arctan(— -^) 

uCv 



(25) 



The plastic modulus is obtained form the modulus of the 
maximal plastic response: 



1 

h 



\da\ 



(26) 



Reciprocally, the plastic response can be expressed in 
terms of these quantities. Let us define the unitary vec- 
tors '0 and . The first one is oriented in the direction 
of and the second one is the 90° rotation of ^/i. The 
plastic strain is written in this basis as: 



(27) 
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FIG. 9. Reduction of the mean coordination number of 
contacts (dotted line). The data have been fitted to a trun- 
cated power law (dashed line). See Eq. (Ea). 
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FIG. 10. Young's modulus. The solid line is the linear ap- 
proximation of E{(£). See Eq. (H). 
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FIG. 11. Poisson's ratio. The dashed line is the quadratic 
approximation of i/(ci). See Eq. (|24[). 
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FIG. 12. Anisotropy parameter. The dashed line is the 
linear approximation of a (d). See Eq. (|^. 
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FIG. 13. Plastic profiles (solid line) and g{e) (dashed 
line). The results for different stress values have been super- 
posed. 

The plastic profiles f{9) and g{9) are shown in Fig. 
|l3| . The first one is approximately the same for all the 
stress states, and can be well fitted to a cosine function, 
centered on the yield direction (/) and truncated to zero 
for the negative values. The last profile depends on the 
stress value, and is difficult to evaluate, because it is of 
the same order as the statistical fiuctuations. However, 
the contribution of g to the total strain response is neg- 
ligible. In order to simplify the description of the plastic 
response, the following approximation is made: 



g{0) « fiO) « [[cos{0 - m = [[0^ 

where [[•]] defines the function 

: X < 
X : X > 



(28) 



(29) 



Now, The fiow rule results from the substitution of 
Eqs. 



and (g8|) into Eq. (jlj): 



J{e)dd 



(30) 



The yield direction and the flow direction have been 
calculated for different stress states. The results are 
shown in Fig. Both angles are quite different, which 
is a clear deviation from Drucker's normality postulate 
p3[ . Indeed, many experimental results on soil deforma- 
tion [Q have confirmed that these angles are completely 
different. Thus Drucker's postulate is not fulfilled in the 



deformation of granular materials, and the main reason 
for that is the rearrangement of contacts on small defor- 
mations which are not taken into account in this theory. 
On the other hand, all the sliding, opening, and other 
micro-mechanical rearrangements can be well handled in 
the discrete element formulation, which results more ad- 
equate to describe the soil deformation. 
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FIG. 14. The How direction and the yield direction of the 
plastic response. Solid lines represent a linear fit. 

The material constants are evaluated from the depen- 
dence of the plastic quantities on the stress: the yield 
direction and the flow direction can be roughly approxi- 
mated by straight lines: 

= 00 + 

P 



V' = -00 +-00- 

p 



(31) 



The four material parameters 0o — 46° ± 0.75°, — 
88.3° ±0.6°, Vo = 78.9° ±0.2° and V'o = 59.1° ±0.4° are 
obtained from the linear fit of the data. On the other 
hand, Fig. ^ shows that the plastic modulus depends 
on the stress through a power law relation: 



h = hn 



1 - 



q 



1 



M Po P 



(32) 



There are four additional material parameters: The plas- 
tic modulus ho ~ 14.5±0.05 at q = 0, the Mohr-Coulomb 
friction coefficient /i* (see Eq. (|l^)), and the exponents 
?7 = 2.7 ± 0.04 and = 0.981 ± 0.002 . 

The plastic limit surface is given by the stress states 
where the plastic deformation becomes infinite. Accord- 
ing to the flow rule ( Eq. (pQ) ), it is found looking for 
the stress values where Eq. (|32D vanishes: 



Po \Po 



(33) 
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It is important to point out that the failure surface 
-given in Eq. (|3|)- does not correspond to the plastic 
limit surface. Actually, This matter has already been 
discussed in the framework of the Hill's condition of in- 
stability the bifurcation analysis [Q, which predict 
that the instability should be reached strictly inside the 
plastic limit surface. 
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FIG. 15. Plastic modulus. The solid line is a power law fit 
with respect to the variable X — 1 ~ {p/vo)* ll {t^* Vo)- 

V. CONCLUDING REMARKS 

The elasto-plastic response of a Voronoi tessellated 
sample of polygons has been calculated in the case of 
monotonic and quasi-static loading. It can be written in 
a simple form as: 

de~=D(d)da + IM!M^ (34) 
h 

The plastic response reflects the non-associated fea- 
tures of realistic soils. Here the yield direction and flow 
direction are linearly related to the ratio q/p, and the 
plastic modulus obeys a power law relation with a weak 
pressure dependence. 

The classical parameters of elasticity - the Young's 
modulus and the Poisson's ratio - are not material con- 
stants, because they depend on the internal damage pa- 
rameter. Therefore, Eq. ( |3^ ) is not complete, and it is 
necessary to include a constitutive equation that relates 
the internal damage to the external load. By focusing on 
the details of the dynamic of the micro-contacts, signif- 
icant progress may be made in the macroscopic descrip- 
tion of the deformation. 

The elasto-plastic response leads to the identification 
of three different regimes which are shown in Fig. 
Zone I corresponds to the isotropic regime, characterized 
by small plastic deformations and a linear elastic regime. 



In the zone II open contacts are detected, which must 
be taken into account in the calculation of the non linear 
elasticity. Zone HI corresponds to unstable states so that 
the stress-strain relation can not be calculated here. The 
extrapolation of the strain response in this region shows 
that the plastic strain must have a finite value just before 
the instability is reached. 

The above observation leads to the open question of 
the nature of the failure . Numerical simulations on 
strain controlled tests show that strain localization is the 
most typical mode of failure. The fact that it appears be- 
fore the sample reaches the plastic limit surface suggests 
that the appearance of the instability is not completely 
determined by the macroscopic state. 

The role of the microstructure on the strain localiza- 
tion has been intensely studied in the last years |po| , pT| . 
Future work is the creation of samples with different 
granular textures -for example, changing the void ratio 
distributions and the polydispersity of the grains-. Then 
we can deal with the question of how a change in the 
microstructure affects the elasto-plastic response and the 
strain localization. 
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